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Instanton methods, in which imaginary-time evolution gives the tunneling rate, have been widely 
used for studying quantum tunneling in various contexts. Nevertheless, how accurate instanton 
^ ^ methods are for the problems of macroscopic quantum tunneling (MQT) still remains unclear be- 

■^1^ , cause of lack of their direct comparison with exact time evolution of the many-body Schrodinger 

' equation. Here, we verify instanton methods applied to coherent MQT. Specifically applying the 

quasi-exact numerical method of time-evolving block decimation to the system of bosons in a ring 
lattice, we directly simulate the real-time quantum dynamics of supercurrents, where a coherent 
oscillation between two macroscopically distinct current states occurs due to MQT. The tunneling 
^ ' rate extracted from the coherent oscillation is compared with that given by the instanton method. 

We show that the error is within 10% when the effective Planck's constant is sufficiently small. We 
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also discuss phase slip dynamics associated with the coherent oscillations. 
PACS numbers: 03.65.Xp,03.75.Kk, 03.75.Lm 



^ ' Tunneling is one of the most fundamental concepts derived from quantum theory and is essential for understanding 
, enormous variety of phenomena in different fields of physics, such as high energy, condensed matter, and atomic 
' physics. The list of such phenomena includes the a-decay of nuclei [ij, tunneling between vacuum states in quantum 
cosmology [11,131 ^'^'^ chromodynamics 0, H, @|i macroscopic quantum tunneling (MQT) in quantum gases 0, [1] and 
■ condensed matter [l^l , and also includes potential applications in quantum information 11 1 . 
C , Instanton methods are general schemes describing quantum tunneling within a semiclassical approximation [Bl. [il. [T^ . 

O . . 

classical equations of motion in imaginary-time coordinate allowing one to obtain an analytical expression for the 
tunneling rate. The instanton methods are closely related to the Langer's formalism of decay of metastable states due 



to thermal fluctuations 13|. Given the versatility and utility of the instanton methods, it is important to examine how 
■ accurately they predict the actual tunneling rate. We note that this question is not entirely trivial. For example, for 
0^ , applicability of the Langer's formalism it is important that the thermal bath (which can be a part of the macroscopic 
' system) is big enough to provide sufficient energy necessary to overcome the barrier separating metastable and stable 
phases. 

00 ■ For single particle problems, especially in one-dimension, the instanton methods can be easily verified because 
\ the corresponding Schrodinger equation can be solved numerically with arbitrary precision. On the other hand, such 
numerical verification of the instanton methods is usually very difficult for complex systems consisting of many degrees 
of freedom, such as MQT and tunneling decay of the false vacuum. Alternatively, in the context of the current-biased 
Josephson junction, where the phase difference between the two superconductors is regarded as a macroscopic quantum 
variable, experiments have been extensively compared with the theory of MQT [l3 Il5|. It has been shown that the 
; ^ ' experiments and the theory are in agreement to the extent that the instanton method provides an estimate of the 
\ order of the magnitude of the tunneling rate. However, this comparison is inevitably limited by the experimental 
uncertainty, which arises from the fact that the theory uses phenomenological parameters extracted from separate 
experiments. For rigorous verification of the instanton methods, it is necessary to make direct comparison of their 
predictions with the first principles many-body simulations in a complex system. 

In this work we study MQT of supercurrents of bosons in a one-dimensional (ID) ring lattice to examine the 
validity of the instanton method applied to MQT. Using the time-evolving block decimation (TEBD) method 1^, 17 1, 



we perform first principles simulations of the real-time dynamics of the corresponding Bose-Hubbard model. In the 
regime where the energies of two macroscopic states with different winding numbers are degenerate, we show that the 
supercurrent exhibits coherent oscillations. These oscillations are accompanied by phase slips which result in sudden 
change of the winding number characterizing the supercurrent. The tunneling rate is accurately extracted from the 
period of oscillations, while it is also calculated by the instanton method in the quantum rotor limit corresponding 
to large filling factors d, [l^] ■ We are thus able to compare the numerical TEBD results with the prediction of the 
instanton method with no ambiguity. Our main finding is that the error of the instanton method is within 10% when 
the effective Planck's constant is sufficiently small. We also find that the coherent oscillations of current persist even 
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FIG. 1: Sketch of the quantum dynamics of supercurrents in the effective potential obtained from the Bose-Hubbard Hamiltonian 
Eq. (m for 9i = n/L (a) and 9i = (b). The blue solid lines sketch the energy landscape versus the current velocity v. The 
black circles represent the quantum state. Note that the plot represents a sketch of an actual process occurring in the multi- 
dimensional phase space. 



between the degenerate states with winding numbers different by two. Such process corresponds to the dynamics 
associated with coherent oscillations of double phase-slips. 



MODEL 



We consider a system of N bosons at zero temperature confined in a homogeneous ID ring lattice of L sites. Recently, 
such a system has been experimentally realized in the context of quantum gases [Toj . We assume a sufficiently deep 
lattice so that the tight-binding approximation is valid. Then, the system is well described by the Bose-Hubbard 
model 

L L 

H = -J^(e-*«6]fe,+i -fh.c.) + - - 1). (1) 

J = l 3 = 1 

where S^+i = hi, reflecting the periodic nature of the ring lattice. The field operator {hj) creates (annihilates) a 
boson on the j-th site, and hj is the number operator. J is the hopping energy and U the onsite interaction. The 



phase twist 9 can be controlled by rotating the lattice [2l|, \2f^ (or equivalently by writing the Hamiltonian in the 
rotating frame). In the case of commensurate fillings, where the filling factor v = N/L is integer, the Bose-Hubbard 
model exhibits a quantum phase transition from a superfluid to a Mott insulator as U / J is increased. Since our 
interest is in the dynamics of supercurrents, we focus only on the superfluid regime throughout this paper. 



SUPERCURRENT DYNAMICS 



For pursuing our main goal of examining the validity of the instanton method, it is imperative to reveal basic 
properties of the quantum dynamics associated with Eq. ([1]). In this section, we confirm that a supercurrent flowing 
through the ring lattice actually exhibits MQT during the real-time evolution. In the next section, we will compare 
the period of oscillations extracted from the MQT dynamics with that obtained by the instanton method. 

To treat the quantum dynamics, we use the quasi-exact numerical method, TEBD 11611 . which is conceptually 
equivalent to the well-known time-dependent density matrix renormalization group 23|, |2J] . This method allows us 
to compute accurately the evolution of many-body wave functions of ID quantum lattice systems. Recently, TEBD 
has been successfully adopted by one of us to a system with periodic boundary conditions 17| . In order to study how 
a supercurrent behaves as a function of time, one needs to prepare a current-carrying state as an initial state of the 
real time evolution. For this purpose, setting 6 = 6q = 27m/ L, we first perform the imaginary time propagation for 
Eq. ID), which provides a current-carrying state with the winding number n. At t = the phase twist is suddenly 
shifted to 9i = 7r(2n — m)/L, m > 1, so that another state with the winding number n — m is exactly degenerate 
with the initial state. We then simulate the dynamics in this system propagating the initial state in real time. First 
we analyze the situation where the filling factor and the initial winding number are equal to unity: ly = l,n = 1 and 
investigate how the time evolution of supercurrents depends on the parameters of the model U/J and di. 

Let us start with the simplest case, 6i = tt/L, where two macroscopically distinct states with winding numbers 1 
and are the degenerate lowest energy states. This situation is analogous to a superconducting flux qubit realized in 
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a superconducting quantum interference device (SQUID), where two flux states with different winding numbers are 
degenerate producing coherent Rabi oscillations [ll|. Likewise in our case we expect coherent oscillations between 
the two degenerate states via MQT as sketched in Fig. [TJa). 

To demonstrate this, we first calculate the time evolution of the current velocity v given by 

- = |^E(^.Vi-i^-)' (2) 

i 

where d is the lattice spacing. When U / J ^ 1, the velocity is almost constant, i.e., the supercurrent is persistent. 
In contrast, when U / J is sufficiently large, e.g. U / J ~ 2.5, quantum fluctuations are strong enough to kick the state 
out from the one of the minima, and the superfluid coherently oscillates between the states with velocities v{t = 0) 
and as shown in Fig.[2l^a). The period of these oscillation decreases monotonically with U / J . 

To confirm that these oscillations are due to quantum tunneling between two macroscopically distinct states, we 
next calculate the overlap |(<I'„|5'(t))p of the wave function with the ground state |<I>„) of the Hamiltonian ([T]) with 
9 = 2-1771/ L, and the momentum occupation n(p, t) — {bj,bp), where bp ~ L^^^^ S^e"*^'^™/^. In Fig. (He), we show 
the overlaps with |$o), and The overlaps and |($o|*(i))P are well approximated by the time 

dependence cos^ijTt/T) and sin^(7rt/T), respectively, where T is the period of oscillations. Hence, the wave function 
is approximated by a macroscopic superposition of the states with n = 1 and ?i = (Schrodinger cat state) as 

\^{t)) c cos + ^ sin (%] |a>o). (3) 



J ' ' \T ^ 

In Fig. [2jb), we show the momentum occupations for p = 1,0,-1, which behave almost identically to the overlaps, 
again justifying validity of the cat state description. We note that the similar cat state dynamics has been found also 
for quantum vortices in anisotropic traps ^25] and supercurrents in two-color optical lattices (26| . 

We next consider the case of 9i = 0, where |$i) and are degenerate. In this case, there arc two possible 

scenarios of the fate of the supercurrent: (i) The supercurrent decays towards the zero momentum state creating 
excitations, (ii) It coherently oscillates between and as sketched in Fig. [T]^b). Previous theoretical work 

on the supercurrent decay anticipated the first scenario to calculate the lifetime of the metastable state using the 
instanton method [28l . |29| . It is very likely that this scenario is indeed realized when the differences in winding 
numbers of Oi and 9q is large. In contrast, it is found in our numerical simulations that the second scenario mainly 
dictates the supercurrent dynamics as seen in Figs. [DJd) and (f). The supercurrent exhibits a coherent oscillation 
between states with velocities v{t = 0) and —v{t = 0) with rapid wiggles. If these wiggles are ignored, then the 
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FIG. 3: (a) Time evolution of the average phase difference 'fi{r,t) for L = iV = 16, U/J = 2.5, and 9i = n/L. The phase jumps 
by 2n at the boarders between the bright and dark regions, (b)-(e) Snap shots of 'fi{r,t) for several values of t. 
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FIG. 4: (a) Time evolution of the average phase difference tf{r,t) for L = N = 16, U/J = 2.5, and (^i = 0. (b)-(e) Snap shots 
of if{r, t) for several values of t. 



wave function is well approximated by superposition of the states |<i>i) and The zero momentum occupancy 

n{p = 0, t) (blue dashed line in Fig. [2]) oscillates in time with the same frequency as the wiggles in the n(jj = ±1, t) 
while the overlap of with |$o) always remains zero. This means that the wiggles come from the coupling with 

the excited states with winding number and that such states contribute to the wave function in addition to 
and 

Since during the coherent oscillations between the two degenerate states, the winding number changes from 1 to 
(or to —1), one expects emergence of the phase slip associated with these oscillations. To reveal the phase slips, we 
calculate the time evolution of the average phase difference between the j-th and (j+r)-th sites, (p{r, t) = a,vg{(p-hj+r))- 
Notice that the phase difference is independent of j because of the homogeneity of the system. In Fig. [3l we show 
ip{r,t) that corresponds to the dynamics depicted in Figs. [2ta)-(c). At t = (Fig. [3^b)), ip{r,t) linearly changes with 
r as (p{r,t) = 2Trr/L corresponding to the winding number n = 1. As time evolves, a phase kink develops around 
r ~ L/2 and it becomes ^ tt at t = T/A (Fig. E^c))- Immediately after t ~ T /A, the phase jumps by 27r and the 
winding number changes to n = as seen in Fig. [2Kd) . 

In Fig.m we show ip{r,t) that corresponds to the situation shown in Fig.[2l^d)-(f), where the supercurrent oscillates 
between the states with n ~ 1 and n = — 1. As t increases, two phase kinks develop; which are localized around 
r = L/A and r = 3L/4. Both phase kinks are tt at t = T/A as shown in Fig. W[c). When t exceeds T/4 (Fig. [Hd)), 
the phase jumps by 27r in the two regions r < L/A and r > ZL/A so that the winding number changes to n = — 1 by 
losing the phase of Air in total. It is worth stressing that this "double phase slip" occurs without passing through a 
state with n = 0. We note that there is no direct connection between the phase slip in real time and that in imaginary 
time [1^ (see also Supplementary Information, Section III) . The dynamics in real time reflects the behavior of the 
average phase difference ip{r,t), which comprises phase slips occurring at different times in different sites. This phase 
slip can be extracted from the superposition of two macroscopically distinct states with different winding numbers. At 
the same time the phase slip in imaginary time develops "instantaneously" during underbarrier tunneling in contrast 
to the phase kink in real time that exhibits a sinusoidal oscillation and develops gradually. Nevertheless the similarity 
between the shape of the phase slip in Figs. Q and ([4]) and the expected shape of the kink in the instanton solution 
is very appealing. 

In the above calculations, we took a relatively small system size L = 16 and unit filling {N = L). Increasing the 
number of sites up to L = 48, we checked that the basic properties of the supercurrent dynamics mentioned above do 
not change. In Fig. [5l we show the period T of the coherent oscillations as a function of L on a log-log scale. There 
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FIG. 5: Period T versus the number of sites L for 6\ — tt/L (a) and 9i = (b). The fiUing factor is fixed to he v = 1. The 
plots are on a log- log scale. 



we clearly see that the period monotonically increases with L following a power law, T cx L". Since the current / at 
a fixed winding number is inversely proportional to L this implies that the frequency of oscillations scales as a power 



of the current. This effect is similar to the situation happening in 2D superconductors at finite temperatures [27[, 
where the supercurrent dissipation rate coming from vortex unbinding also scales as a power of the current. We 
also note that the commensurability of the filling factor is crucial for the coherent supercurrent dynamics. Only in 
case of commensurate fillings, the two states |$i) and |<i>o) (or are coupled through the Umklapp scattering 

process [2l|, |22| and the coherent oscillations can occur. 



COMPARISON WITH THE INSTANTON METHOD 



Having established that the supercurrent dynamics exhibits coherent oscillations caused by MQT, we now compare 
the frequency of oscillations calculated by the instanton method with the TEBD results. For this purpose, we choose 
the situation of 0i = tt/L, characterized by a single phase slip dynamics. To seek a simple analytical expression of the 
energy splitting between the two current states A, which is the same as the oscillation frequency, we assume Uv^ J 
and V ^ \. In this case the Bose-Hubbard model can be mapped onto the 0(2) quantum rotor model [1, This 
limit also describes a regular array of coupled Josephson junctions. We have confirmed that the quantum rotor model 
indeed gives a very accurate value of A when the filling is large: v > 1000 (see Supplementary Information, Section 
II). In the quantum rotor model the phase space variables are the superfluid phases on each site and the conjugate 
momenta, corresponding to the fluctuations of the number of particles. So overall the phase space consists of 2L 
variables. Applying the instanton method to the c^uantum rotor model, we obtain the energy splitting expressed as 



# = W:^expf-fi), (4) 



where = V vJU is the Josephson plasma energy and /ic = ^/UjiyJ) is the effective Planck's constant. Note that 
he <^ \ deep in the superfluid regime where number and phase can be approximately treated as classical variables. At 
/le ~ 1 the quantum fluctuations become important and can even drive the system to a different insulating phase [ist . 
Since the instanton method is a generalization of the WKB semiclassical approximation 0, , the expression for the 
energy splitting ^ is supposed to be accurate when h^/sj is sufficiently small. We note that the quantum rotor model 
has a clear advantage over original Bose-Hubbard model that the instanton action sj and the coefficient K do not 
depend on U/ J and u (see Methods and Supplementary Information for specific expressions of s/ and K). In other 
words, the dependence of A/Ej on U/J and v comes only through a single parameter h^. Using explicit calculation 
(sec Methods) we obtain sj = 7.363 and K = 3.06 if we set the number of sites to be L = 8. 

We can also extract the energy splitting A from TEBD simulations by fitting the the overlap f{t) = |((f>i|\I'(i))p 
(like in Fig.[2l^c)) using the function 

f{t)=Bcos^f^t\+C, (5) 



where A, B, and C are the free parameters. In Fig. [Blja), we show the energy splitting A versus he calculated by the 
instanton method (blue solid line), and by TEBD for the filling factors v = 1000 (red squares) and ly = 10 (black 
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FIG. 6: (a) Energy splitting A/Ej as a function of the effective Planck's constant hs = sJlJ I {yJ) for L = 8. The blue solid 
line represents the result by the instanton method corresponding to Eq. Q with si = 7.363 and K = 3.06. The red squares 
and the black circles are the TEBD results for v — 1000 and u = 10, respectively, (b) Ratio |Atebd — Ains|/ATEBD of the 
difference between the TEBD and instanton results as a function of ha- 

circles). It is evident that for ly = 1000 and ft-c sufficiently small the instanton and TEBD results agree very well. 
To quantify the error of the instanton method, in Fig. [D^b) we show the relative difference between the two results: 
|Atebd — Ainsl/AxEBD- For i/ = 1000 (red squares), as decreases, the error also decreases such that it is within 
10% when < 0.7. It is hard to push the calculation to even smaller values of because of exponential sensitivity of 
the period of oscillations to the effective Planck's constant. Nevertheless our results allow us to make the conclusion 
that the instanton method can provide quantitatively accurate prediction for the tunneling probability when he/sj 
is sufficiently small. At the same time, the error for = 10 is significantly larger than that for v — 1000. Moreover 
the error does not even monotonically depend on h^. This clearly means that at this filling the quantum rotor model 
gives only qualitative description of the tunneling process. 

SUMMARY 

We analyzed quantum dynamics of supercurrents of one-dimensional lattice bosons in a ring. In particular, our 
focus was on the coherent oscillations between the two degenerate current states via macroscopic quantum tunneling 
(MQT). The period of these oscillations T is related to the energy splitting A induced by the tunneling as T = 27r?i/A. 
We calculated A both simulating real-time dynamics using the time-evolving block decimation (TEBD) method and 
within the imaginary time instanton method. We showed that the result of instanton calculation is in very good 
quantitative agreement with the TEBD result when the effective Planck's constant is sufficiently small. This 
agreement verifies the instanton method applied to coherent MQT involving many collective variables. 

We also want to emphasize that the success in applying TEBD (or equivalently the time-dependent density matrix 
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renormalization group) to MQT problems opens up new possibilities to analyzing macroscopic tunneling phenomena. 
In particular, (i) TEBD allows us to precisely calculate the energy splitting even for large h^, where the instanton 
method fails, (ii) TEBD provides time evolution of the many-body wave function, from which one can calculate various 
quantities, for example different correlation functions. The first advantage (i) is crucial for quantitative simulation 
of experiments (e.g. in cold gases), where it is easier to work in the regime of larger h^. and shorter periods to 
avoid various effects of decoherence like particle losses. Moreover, the second advantage allows one to reveal detailed 
processes of MQT in real time. As an example, we have analyzed the time evolution of the phase-phase correlation 
functions and revealed the existence of the phase slips associated with the coherent oscillations, which can be detected 
in experiments. One can extend this analysis to study higher order correlation functions to e.g. detect shot noise of 
phase slips or even their full counting statistics (30| . 

I. D. thanks M. Nishida, A. Nunnenkamp, T. Nikimi, S. Kurihara, and Y. Kato for valuable comments and dis- 
cussions. I. D. acknowledges support from a Grant-in-Aid from JSPS. I.D. is grateful to Boston University visitors 
program for hospitality. A. P. was supported by AFOSR YIP and Sloan Foundation. 



METHODS 



Here, we outline instanton derivation of the energy splitting ([4]) and give specific expressions of s/ and K. For 
additional details of the derivation of these expressions we refer the reader to the Supplementary Information, Section 
III. In the limit when Uv ^ J and 3> 1, number fluctuations are significantly suppressed and the Bose-Habbard 
model Eq. ([T|) can be mapped onto the 0{2) quantum rotor model [1], described by the effective action 



df 



ld(j) d(j) 
2df ' df 



where (j) is an L-dimensional vector defined as 
and the potential is 



(6) 



(7) 



v{$) = 5^y,(0,+i,0,) 

L 

= ^ -2 cos - 4>j ~ 6) . 



(8) 



is the phase of particles at the j-th site, and f = tEj/H and /3 — Ej/{k^T) denote the imaginary time and the 



inverse temperature in the Josephson plasma energy unit. We remind that Ej = \J vJU. 

Extremizing the action, by setting 8s — 0, we obtain the classical equations of motion for the phases: 



9f2 



= 2 sin — — 6) — 2 sin ((^j — (jij-i — 



(9) 



In order to calculate s/, we numerically find the instanton solution 0(f) = </>/('?) of Eq. which connects the two 
degenerate states with different winding numbers: 



0,(-;9/2) = ^ - ^ ("l + i ) , 0^.(^/2) = 0. 



(10) 



(see Fig. 7 of the Supplementary Information). Once this solution is obtained, we obtain s/ by substituting 0(f) 
$i{f) into Eq. ([6]): 
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For L = 8, Eq. ^ gives s/ = 7.363. 

In turn the prefactor K is given by [1, Q 



K = 



n„.A 



(0) \l/2 



where Am's and Am'''s are the solutions of the following eigenvalue equations: 
and 



Here the L-dimensional vectors 



and 



6n = (Cl,m(T), . . . ,^j,m(T), . . . ,CL,m(f)) , 



obey the orthonormalization condition 

The L X L matrices yVJ and jCi are defined as 



and 



where = S^^. V, 



■^SS - hk 



(12) 



(13) 



(14) 



(15) 



(16) 



(17) 



(18) 



(19) 



^ ^. Notice that (L + l)-th and 0th sites are equivalent to 1st and L-th sites, respectively, reflecting 



4=0 



the periodic boundary condition. For L = 8, Eq. (|12p gives K = 3.06. 
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SUPPLEMENTARY INFORMATION 



I. TEBD FOR LARGE FILLING FACTORS 



In this section, we present an idea of adopting the time-evolving block decimation (TEBD) method to the Bose- 
Hubbard model when the average number of particles per site v (or the filling factor) is large. The key of the idea is 
that in addition to the upper bound, the lower bound for the occupation number of particle per site is introduced in 
order to significantly reduce the size of the local Hilbcrt space. This idea is crucial because the quantitative comparison 
of the TEBD results with the results of the instanton method based on the quantum rotor model is possible only for 
very large v > 1000 (see Sec. II). 

Let us consider a system described by the one-dimensional (ID) Bosc- Hubbard model with L lattice sites. Spanning 
the Hilbert space of the whole system by a product of local Hibcrt spaces of dimension d, a many-body wave function 
of the system is expressed as 
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FIG. 7: Occupation probability P{n) (in the log-scale) of the local Fock state \n) in the ground state of the untwisted Bose- 
Hubbard model with L = 8, = 10, and U/J = b. Here L is the system size, u is the filling factor, U is the onsite interaction 
and J is the hopping energy (see Eq. (1) in the main text). 
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In the TEBD algorithm [16|, coefBcients j2_...jj^ are decomposed in a particular matrix product form as 

X 

•-Jl J2,--- Jl, — Qi ^Qi-^ QiQ2^Ct2 CtL-2 QI,-2QI,-l^CtL-l CtL-1' V 

ai ,...,Q!t-l = l 

The vector Aq| represents the coefficients of the Schmidt decomposition of j^I') with respect to the bipartite splitting 
of the system into [1, . . . , Z — 1, ^] : [Z + 1, Z + 2, . . . , L]. The tensors F's constitute the Schmidt vectors together with 
the A-vectors. x is the number of basis states, which is taken to be sufficiently large so that the error due to this 
truncation is nearly equal to zero. In the typical calculations in the main text, it ranges from x = 100 to x = 200. 

Usually dimension of the local Hilbert space corresponding to a single site is chosen as d = nmax + 1, where n^ax 
is the maximum number of particles per site. It is spanned by the basis set, {\n = 0), |1), . . . , In^ax — 1), |"-max)}- 
While in principle nmax is equal to the total number of particles in the system, taking much smaller rimax provides 
converged results in practice. For instance, for accurate determination of the zero temperature phase diagram of the 
Bose Hubbard model at unit-filling, Umax = 5 (d = 6) is sufficient [3l|. At large filling factors, however, this choice 
of the local Hilbert space basis makes computations extremely expensive, because the computational cost in TEBD 
scales as Ld^x^ fl6| . To solve this problem, in addition to rimax, we introduce the minimum number of particles per 
site riinin and span the local Hilbert space by the basis set, {\n = nmin), I'T-min + 1), • ■ ■ , I'^-max — 1), I'T-max)}, and thus 
d — ?T.max ^ ?^min + 1- In thc parameter region of U /{vJ) ~ 1, where our TEBD simulations are carried out, setting 
'T-max = + 5 and riniin = 1^ — 5 corresponding to d = 11 is sufficient for the convergence regardless of the value of v. 
For instance, in Fig. [71 we plot the occupation probability P(n) of the local Fock state \n) for L = 8, = 10, and 
U / J = 5. It is evident that P{n) exponentially decays as n deviates from its average v = \Q and that P{n) for n > 15 
and n < 6 is less than 10"^. This justifies this truncation scheme for practical calculations. 



II. EFFECTIVE ACTION FOR THE PHASE SLIP PROBLEM 



In this section, we explain the mapping of the Bose-Hubbard model onto the 0(2)-quantum rotor model. For this 
purpose, we start with the grand canonical partition function, 

S[h*,h] 



Z = / Vh*Vhcy.^l - 



d 



(22) 



where the action S[b*, b] for the Bose-Hubbard model (see Eq. (1) in thc main text) is given by 

L 

/ 2 

dT 



b*{r)h—b,{r) - J {b*{T)b,+,{T)e-'^o + b*^,{T)b,{T)e'') 



U 



-b*{T)b*{T)b,{T)b,{T) - M6*(T)6,(r) 



(23) 



where U is the onsite interaction, J is the hopping energy, v is the filling factor, fi fa Uv is the chemical potential, 
and 9 is the phase twist. For convenience we introduce finite small temperature T corresponding to the inverse 



temperature (3 = {UbT)^ . In thc end of calculations we will take the limit T ^ 0. Inserting bj 
is rewritten as 

L „M 
S'[n, 4>\ =y I dT hn-j 

j=l-'-— 



, thc action 



d(j)j 1 drij \ 

a7 + 2^^ 97 ) ~ 2VwrJcos(0,+i 



(24) 



Splitting thc number of particles per site into its average and fluctuation as nj ~ v + Snj , assuming that i> is integer 
and that Uv ^ J and v ^ Srij , we find that the action is then approximated as 

J, hfi 

d(j)j 



= V / ' dT 

3 = 1-^ ~ — 



ihSn 



2iyJ cos {4>j+i 



2 ^ 



(25) 



Since Eq. (pS)) contains only the linear and quadratic terms with respect to number fiuctuations Snj , these degrees of 
freedom can be integrated out. Then, the action is described in terms of the phases as 



= E / dT 



2U [~d^ 



— 2iyJ cos ( 



- ^3 ' 0) 



(26) 
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It is convenient to change variables 



so that Eq. (f26|) is rewritten as 



where s is the dimensionless action 



M = ^ / dr 

j-^l J -13/2 



2 cos( 



(27) 



(28) 



(29) 



and (3 — /3V vJU . From Eqs. ((22)) and p8|) we clearly see that = y/U /{vJ) plays the role of the effective 
dimensionless Planck's constant for this problem with /lo corresponding to the classical (Bogoliubov) limit and 
/lo > 1 corresponding to the regime of strong quantum fluctuations. 

Extremizing the action by imposing 5s ~ 0, we obtain the classical equations of motion for the phases (f>j, 



= 2 sin {(pj+i - - 6) - 2 sin (0^ - 4>j-i - 



There are two types of stationary solution of Eq. ([30]) . One is 

2?™ 



L 



(i-1), 



(30) 



(31) 



which describes the current carrying states with the winding- number n. The other is a saddle-point solution with a 
phase kink separating (meta)stable states with different winding numbers: 



+ - 1), 



(32) 



where 



and 



L-1 + 2n 
L-2 



29- 



L-2 



mod 27r, 



27:71 — a 
L-1 ' 



(33) 



(34) 



Notice that in Eq. (|33p the phase kink is assumed to be located at the link between the 1st and L-th sites. The 
magnitude of this kink a is defined within the interval [— 27r,0]. In the limit of the large number of sites L 3> 1 the 
expression for a simplifies: 



-TT 1 



2n 
T 



26 mod 27r. 



(35) 



In particular, in the case n — Q and 6 ~ tt/L. which wc arc interested in, a — — 7r(l — 1/^) ~ — tt and (p — n/L. 

For 9 ~ tt/L, the two current states with windings n = I and n = are degenerate. Quantum tunneling couples 
them and breaks the degeneracy, leading to the energy splitting A between the ground (bonding) state and the first- 
excited (anti-bonding) state. This tunneling process is associated with generation of a "phase slip" or cquivalcntly 
a phase kink. In imaginary time evolution the virtual kink forms during the imaginary time evolution of the phase 
between the two current states. If the state of Eq. ()3ip with n = 1 is prepared initially, the many-body wave function 
coherently oscillates with the period 27rft/A between the states with n = 1 and n = 0. It is well known [I4I that the 
energy splitting can be expressed as 



(36) 
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FIG. 8: Ratio of the energy splitting A to the Josephson energy Ej — yvJU as a function of v for U / i^vj) = /i^ = 1 (a-) a-nd 
U/{yJ) = 0.5 (b). The solid lines represent the value of A/i5j at = 5000. The dashed lines are the prediction of the instanton 
method. Notice that both U and v are changed for a fixed value of J such that remains constant. 



where 



and 



A- 



Zi = 



(0) 



(37) 



(38) 



Notice that ^^^^ Vcj) denotes the path integral over trajectories containing a single instanton, while J^^^ Vcj) is the path 
integral containing zero instantons. According to the instanton method 0, the energy splitting of Eq. ([55]) is well 
approximated by 



A ~ 2LKEj 



■ exp 



(39) 



where s/ denotes the action for the instanton solution. K is the constant we define below, and i?j = \/ vJU is the 
Josephson plasma energy. It is worth stressing that s/ and K do not depend on v and U / J, but depend only on L so 
that A/Ej depends on U, J, and only through he- In the following section we will present a derivation of Eq. ([55]) 
first evaluating it approximately and then exactly and will give the explicit form of the coefficient K . 

As mentioned above, the mapping onto the quantum rotor model is justified when Uv ^ J and 3> 1. For 
unambiguous comparison between the TEBD and instanton results, we need to specify quantitatively the parameter 
region where the quantum rotor model is valid for the calculations of the energy splitting. It is clear in Eq. ([M)) as 
V increases at a fixed value of /le, A/ Ej should saturate at a constant corresponding to the quantum rotor limit. As 
we show in Fig. [51 where we plot A/Ej versus for two different values of h^, this is indeed the case. Note that this 
ratio A/Ej becomes independent on i/ only at very large fiUiwng factors i/ > 1000. We also want to point that for the 
smaller value of he the larger the filling factor is required for the convergence. Since our quantitative analysis in the 
main text is focused on the region of U/{uJ) > 0.5, the quantum rotor model is sufficiently accurate for v — 1000, 
which we use in practice. 
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III. INSTANTON METHOD FOR THE QUANTUM ROTOR MODEL 



A. One collective variable 



Since the instanton method in the presence of many degrees of freedom is quite compHcated, we first use a simpler 
model, which is reduced from the quantum rotor model by assuming that the phase slip is described by only a single 
collective variable. This simple model, which represents a variational estimate of the full result, is useful to understand 
basic ideas of the calculation. Later we will generalize the result to the phase slip described by two degrees of freedom 
and finally show the complete instanton solution of the full problem. In the regime of validity of the quantum rotor 
model, the healing length ^ = d\/23 j (i^U) is much shorter than the lattice spacing d. Hence the phase slip that 
develops during the tunneling process should be localized with in a few sites. Without loss of generality, we can 
assume that the phase kink develops at the link between the 1st and L-th sites. In the first approximation we assume 
that the phases along the instanton trajectory satisfy the ansatz, 

cl^,{f) = ^+'p{r){j~l), (40) 

where a denotes the phase difference between the 1st and L-th sites. Its time dependence is found from extremizing 
the effective action (see below) . The remaining phases on other sites are chosen as a simple linear function of the site 
index j with (p = —a/{L — 1) such that the boundary condition (j>L = —a/2 is fulfilled. Substituting Eq. (|40l) into 
Eq. ([29]) we find that the effective dimensionless action describing the system becomes 



= / df 



(41) 



This is nothing but the classical action of a particle with the effective mass M, which depends on the system size 
according to 

M=^^, (42) 
12(i-l)' ^ ' 

moving in the effective potential — Vi(q:), where 

Vi{a) = -2 cos (a - 61) - 2(L - 1) cos(^ - 61). (43) 

The shape of Vi{a) — Vi(0) for 9 ^ tt/L and L = 8 is depicted in Fig. [9l Note that Vi{a) has two global minima 
a = ai = — 27r(l — l/L) and a = a/ = corresponding to the current-carrying states with winding numbers n ~ 1 
and n = respectively. These two minima are separated by a local maximum, a ^ as = — 7r(l — l/i) describing the 
saddle-point solution of Eq. ([55]) . Thus, introducing the collective variable a, the phase slip problem is equivalent 
to tunneling of a single particle in a one-dimensional symmetric double-well potential. The corresponding classical 




FIG. 9: Effective potential Vi{ol) for Q = t^IL and i = 8. 
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FIG. 10: Instanton solution ai{f) and dfai{f) 



equation of motion describing the particle motion in the (inverted) potential —Vi{a) is 



dc 



(44) 



The instanton solution of this equation a{f) ~ aj{f) is the one satisfying the boundary conditions q;(— /3/2) = and 
a{f3/2) = Oif. Such a solution (shown in Fig. llOp contains a kink in the phase a. The instanton solution defines the 
classical trajectory in the path integral of Zi. There is another trivial solution of Eq. (j44p . a{f) = ai (or equivalently 
a{f) = a/), which is the classical trajectory corresponding to the path integral of Zq. 
To calculate the ratio A in Eq. ((37|) (see Ref. Q for more details), we substitute 



a(f) = ai 



(45) 



into Zi and 



a(f ) = ai + 



(46) 



into Zq, where ^m's and fm^'s are complete sets of real orthonormal functions obeying the following eigenvalue 
equations: 



1 d^Vi 



92 



9f2 M 



^m(T) = A„^„(f) 



and 



(47) 



(48) 



with = M ^d'^Vi\a=ai- Neglecting the terms higher than the second order with respect to yHi^JM, A is approxi- 
mated as 



exp - — 



rsA !■■■! n,j2^)-'/'dc™ exp [-1 A™c^] 



'° ^ / • ■ ■ / n™(2^)-^/2rfc„ exp [-IE™ a5;^^c2 
where s/ is the action of the instanton solution given by 



(49) 



si = dfM 



( 

\ df 



(50) 



To carry out the integrals with respect to Cm's in Eq (|49p . it is important that due to the translation invariance of 
the instanton solution in the imaginary time, Eq. (|47p possesses one solution ^0 with the eigenvalue Aq = 0. For this 
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zero mode, the integral in Eq. is formally divergent. To solve this problem, one needs simply replace J dco with 
^/sj/hc / dr leading to 

In the above discussion, we assumed that the phase kink develops at the link between the 1st and i-tli sites. In total 
there are L independent possibilities for the kink. Note that because we are dealing with a discrete system, there 
is no continuous symmetry associated with this degeneracy and thus no additional zero eigenvalue in Eq. (j5ip . All 
instanton solutions centered around different links give identical contribution to Zi . It is therefore only necessary to 
multiply Ahy L before substituting it into Eq. (|36p . Thus, we obtain Eq. ((39)) with the coefhcient 

C(o) \ ^/"^ 
. (52) 
1 hn^f) / 

Now both s/ and K can be straightforwardly found numerically. For the situation of a single collective variable 
described here they are explicitly given in the first row of Table H] 



B. Two collective variables 



After considering a toy single-variable approximation to the instanton solution, in this section, we make the next 
step by increasing the total number of the collective variables to two. Specifically to describe the instanton action 
we take two variables a and (3 describing the phase slip as independent and for the rest use the linear interpolating 
function. The phases of such instanton solution (again centered between 1st and L-th sites) arc described as 

r a(f )/2, for j = 1 

cj),{T) = <^ a(f)/2 + /3(t) + v{f){j - 2), for 2 < j < L - 1 , (53) 
I -a(f )/2, for 3 = L 

where a denotes the phase difference between the 1st and L-ih sites, /3 is the phase difference between the 2-nd and 
1-st (as well as L-th. and {L — l)-th) sites, and ip = —(a -I- 2(3)/{L — 3) is chosen such that the boundary condition 
~ —02 is fulfilled. Substituting Eq. ([53]) into Eq. ((29|) we find that the action is described by the two variables 
a and /3 as 




solution. 
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FIG. 12: Instanton solution xi(f) and yi{f) for 9 = n/L and L = 8. 



where 



_ £^ + 3L + 16 ^ _ (£-l)(£-2) 

^ " " 12(L-3) ' ~ - 3(L-3) ' 

and V2(a,/?) is the following effective potential 

1/2 (a, /3) = -2cos(a-6l)-4cos(/3-6')-2(L-3)cos(v3- 



(55) 



(56) 



It is convenient to perform a linear transformation {x,yY = X(a,/3)*, where X is an orthogonal 2x2 matrix, to 
diagonalize the kinetic energy part of the action leading to 



s[x, y] 



df 



dx 



dy 



(57) 



The shape of V2{x, y) — V2(0, 0) for 9 = tt/ L and L = 8 is depicted in Fig.[TTJ In the potential there are two minima 
corresponding to the current states with n = and n = 1. The classical equations of motion corresponding to this 
action are 



d'^x 
d'^y 



dx 
dV2 



dr"^ dy 



0, 



0. 



(58) 



As in the case of the single collective variable, the instanton solution describes the classical trajectory in the inverted 
potential, which starts from one of the maxima of —V2{x,y) (or equivalently minima of V2{x,y)) at — f ~ 13/2 and 
reaches the other maximum at f = 13/2 through a saddle point as shown in Fig. 1121 The trajectory of the instanton 
solution is indicated by the solid line in Fig. [TlJ Inserting the instanton solution into Eq. ([57)) . we obtain s/. The 
derivation of the coefficient K in Eq. ((5^ is almost the same as that for the single collective variable and we skip it to 
avoid redundancy. The values of both 5/ and K in this two- variable case can be found in the second line of Table [H 
In a similar way one can keep on the number of independent degrease of freedom in the instanton solution. 



C. All degrees of freedom 



As a final step we will explicitly show generalization of the instanton method to the action of the quantum rotor 
model Eq. ([29|) where all phases are treated as independent variable. We will show that the energy splitting is given 
by Eq. (|5ip where the eigenvalue equations (|47p and (j48p are appropriately generalized. For convenience, we rewrite 
Eq. dig) as 
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where 6 is an i-dinicnsional vector defined as 



and the potential is 



L L 



(60) 



(61) 



The classical equations of motion Eq. (j30p have an instanton solution 4>{f) = (j}i{f) that connects two current states 
through the saddle point having a phase kink. We obtain such a solution by numerically solving Eq. pop imposing 
the boundary conditions: 



0,(-/3/2) = ^ - ^ A + i ) , 0^.(/3/2) = 0. 



(62) 



The corresponding instanton solution for L = 8 is depicted in Fig. 1131 Notice that apart from the kink between 
4-th and 5-th sites the remaining phases approximately linearly depend on the site index justifying the single- variable 
variational ansatz made in Sec. . However, because of high sensitivity of the splitting A to especially the value of s 
such ansatz can not be used for accurate quantitative calculations. We intentionally shifted the position of the kink in 
Fig. [13] to the middle of the system for better graphical presentation. For computational purposes it is convenient to 
assume that the link develops between 1st and L-th sites as we did in earlier calculations. Substituting 4>{t) = 4>i{f) 
into Eq. (|59|1. we obtain the instanton action s/. 
To calculate A of Eq. ((37|) , we substitute 

$(f) = ^/(f) + Vfe^y^Cmgm(f), (63) 



into Zi and 



(f) = ^(-/3/2) + c™el°nr), 



(64) 



into Zf), where 



iin = (^1,^(7"), ■ • ■ ,£,j,m{T), ■ ■ ■ ,iL,m{T)) , 



(65) 



(66) 




FIG. 13: (a) Instanton solution 4'i{t) for 6 = n / L and L = 8. (b) Snap shots of 4'i{t) for f = —6 (black diamonds), (blue 
squares), and 6 (red circles). 
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Number of collective 
valiables: m 


Instanton 
action: s/ 


Coefficient: 
K 


1 


7.749 


3.71 


2 


7.396 


4.89 


3 


7.364 


4.41 


4 


7.363 


3.64 


8 


7.363 


3.06 



TABLE I: si and K for several values of the number of collective variables, where L 



The L dimensional vectors |^i's and |^'''s obey the eigenvalue equations 



-M(°^e(r)=A(::'a(r), 



and the orthonormalization conditions 



dr ■ U = 5ur.. j df ■ 4°) = 5i^^. 
The L X L dimensional matrices tW and tW^"^ are determined by the matrix elements 



M 



5j,k 



Q^2 



50? 



d(j)jd(j) 



'j+i 



d(l)jd(j) 



(67) 
(68) 

(69) 

(70) 
(71) 



where = 9? ^ ^. Notice that (L + l)-th and 0th sites arc equivalent to 1st and L-th sites, respectively, reflecting 

^ 0—0 

the periodicity of the system. Neglecting the terms higher than the second order with respect to ^/h^, A is again 
approximated as Eq. (|49p . The derivation of Eq. p9|) with the coefficient K given by Eq. ((52l) is exactly the same as 
the case of the single collective variable and will not be repeated here. The only difference with the single variable 
case is that A^'s and Am^'s are now given by the eigenvalues of Eqs. ([67|) and ([68|l . 



D. Comparison with the TEBD results 



In the main text we compared the energy splitting calculated by the instanton method with all possible degrees of 
freedom (as described in Sec. ) with the TEBD results. It is also instructive to learn how the instanton method is 
improved as the number of collective variables increases. For this purpose, let us now present the comparison between 
the TEBD and the approximate instanton results where only m < L collective variables are treated as independent. 
We take a relatively small system size L = 8. In Table U we first show the instanton action sj and the coefficient K 
for several values of m. Both sj and K approache the exact values corresponding to m = 8 as m increases. We note 
that the action sj for m = 4 is exactly the same as that for m = 8 because the instanton solution is anti-symmetric 
with respect to j ^ L — j, i.e. with respect to the link at which the phase kink develops (see Fig. [13]). In contrast, 
K for m = 4, where the fluctuations are also forced to obey the same symmetry as well, is significantly different from 
if for m = 8. Thus, it is crucial to include all possible fluctuations in order to obtain the correct value of K. 

In Fig. [Mia), we plot the energy splitting calculated by the instanton method as a function of he together with 
that by TEBD for ly = 1000. At first glance, it seems that the results by the instanton method with a single collective 
variable agrees very well with the TEBD results. However, this seeming agreement is rather coincidental as shown 
in Fig. [Ml where we plot the ratio |Atebd — Ains|/Ains of the difference between the energy splittings by TEBD, 
AtebDj and the instanton method, Ai„s. There we clearly see that the error for m = 1 (black triangles) does not 
monotonically decrease with h^, contradicting the basic fact that the instantons should be more accurate at smaller 
he- Except for this case with m = 1, the error decreases monotonically as the number of independent phases m 
increases and the effective Planck's constant decreases. 
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FIG. 14: (a) Energy splitting A/Ej versus the effective Planck's constant he = \/U/{i'J). The dotted, dashed, and solid lines 
represent the results by the instanton method for m = 1, 2, and 8. The dots are the TEBD results for u = 1000. (b) Ratio 
|Atebd — Ains|/Ains of the difference between the TEBD and instanton results as a function of he. 



